Ripple and kink dynamics 
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We propose a relevant modification of tlie Nisfiimori-Oudii model [Phys. Rev. Lett. 71, 197 (1993)] 
for granular landscape erosion. We explicitly introduce a new parameter: the angle of repose Or, 
and a new process: avalanches. We show that the 6r parameter leads to an asymmetry of the 
ripples, as observed in natural patterns. The temporal evolution of the maximum ripple height 
hmax is limited and not linear, according to recent observations. The ripple symmetry and the kink 
dynamics are studied and discussed. 
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I. INTRODUCTION 

The formation of ripples due to the wind blowing 
across a sand bed [1] is a complex phenomenon involv- 
ing various mechanisms such as air turbulence, non-linear 
fluxes of sand, saltation, grain ejection threshold, creep- 
ing and avalanches. 

Different models [2-5] for ripple formation have been 
proposed. Some of those models consider both hopping 
and rolling of grains, and they are able to reproduce trav- 
elling ripple structures. Simulations [6] have also consid- 
ered various additional effects like the screening of crests, 
the grain reptation and the existence of a grain ejection 
threshold. 

The Nishimori-Ouchi (NO) model [3] is able to repro- 
duce a wide variety of different eolian structures: ripples, 
star dunes, barkhanic dunes, etc... The main advantage 
of such a numerical model is that all parameters can 
be easily tuned and the dynamics can be deeply investi- 
gated. In a recent work [7], we have shown that the NO 
model reproduces the complex labyrinthic patterns cre- 
ated in some experiments [8]. However, the NO model 
leads to unrealistic features. Among others, a cube of 
sand can be stable within the NO model. It is possible 
to create aeolian structures with an infinite height. In 
fact, the existence of a critical angle is lacking in the NO 
model. 

In the next section, we will propose a modification of 
the NO model in order to impose a more realistic con- 
straint to granular surfaces. In section 3, we will study 
how the ripple formation depends on the various param- 
eters of the model. In section 4, the dynamics of the new 
model will be studied and will be confronted to various 
situations. Finally, a summary of our findings will be 
given in section 5. 



II. MODELS 

In this section, we present two 2-dimensional models 
of granular landscape evolution. 



A. Nishimori-Ouchi (NO) 

In the Nishimori-Ouchi model [3], two modes of gran- 
ular transport are considered: (i) saltation and (ii) rep- 
tation. The saltation process acts along the x axis. It 
takes grains from a site {x — i, y) and add these grains 
to a site {x,y). The value of the saltation length £{x,y) 
depends on the local slope of the landscape at {x,y). In 
this model, the saltation grains are captured by the sur- 
face when they touch it. The granular flux Qsait due to 
the saltation process can be written as 



Qsait^A N{x',y)dx' , 



(1) 



where Nix, y) is the amount of sand per time unit and 
surface unit at the position {x,y). A is a dimensionless 
constant parameter. 

The NO model also considers a 2-dimensional diffusion 
process, i.e. the reptation. The amount of sand displaced 
at (a;, y) is supposed to be proportional to the local curve 
of the surface at that position. The flux Qrept due to the 
reptation reads 



Qrept — D^h, 



(2) 



with a kind of diffusion coefficient D. The mass conser- 
vation imposes 

dh 



dt 



= -VQ, 



(3) 



where Q = Qsait + Qrept- Therefore, the temporal evo- 
lution of the height h of sand at the position (x, y) reads 



dh{x,y,t) 
dt 



A 



N{x~£,y) 1 



diix,y) 
dx 



N{x,y) 



DAh. 



(4) 
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As in the NO model [3], we have assumed that the 
local slope (along the x direction) mainly controls the 
sand flux Qsait- Experimentally [1], one observes that 
the wind is stronger on the crests. As a consequence, 
the amount of displaced sand is larger around the crests. 
This cfl^ect can be included in the NO model by assuming 



dQsal. 

dx 



5\l + tanh 



dh 
dx 



1 



tanh 



dh 
dx 



(5) 



where the parameters 5 and e give respectively the scale 
factor of the saltation flux and the non-zero asymptotic 
quantity of sand displaced when dh/dx > 0. The salta- 
tion length £ is represented by the factor (tanh(H) -f 1). 
This assumes that i is maximum on the crest and that 
small saltation jumps take place on the downstream part 
of the ripple. 

We have numerically solved the NO model on two- 
dimensionnal square lattices starting from an initially 
plane landscape. At each time step, a randomly cho- 
sen site is submitted to both saltation and reptation pro- 
cesses (Eqs.(4) and (5)). Figure 1 shows the maximum 
height hjnax reached by the proflle as a function of time 
t. Time is expressed in numerical steps (iterations). 
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1. Evolution of the maximum height hmax of the sur- 
face as a function of time t in the case of the NO model. The 
simulation parameters are : S — 25, D = 0.4 and e = 0.3. 
The lattice size is 101 x 101. 



B. Saltation-Creep-Avalanche (SCA) 

Granular materials are characterized by the angle of re- 
pose Or [9]. If the slope of a granular pile exceeds tan(0r), 
an avalanche is created and stabilizes the surface. The 
angle of repose depends on the grain characteristics: size 
distribution, shape, rugosity, and density. Typically, this 
angle 9r is about 37° for sand [1,9]. Since it is an intrinsic 
property of the granular material, the relevant parameter 
Or should be considered in numerical models. 

Let us consider a square lattice with discrete spatial 
variables x and y. The stability condition with respect 
to the angle of repose Or reads 



^2^+1 1 < tan(^^J.) 



\hy — hy+i\ < tan( 



(6) 



where the left hand side is the discrete counterpart of the 
spatial derivative of the height of sand h at the position 
{x,y). If Eq.(6) is not verifled for a speciflc direction {x 
or y), an avalanche occurs as illustrated in Figure 2. In 
this case, the excess of sand E = tan(6') — tan{Or) (with 
tan(0) — with i = x or i = y depdending on which 
condition (6) is not verifled) is supposed to be equitably 
distributed between x and x + 1 ii the avalnche occure 
in that direction. The same applies for the y direction. 
Taking those considerations into account, Eq.(4) becomes 



dh{x,y,t) 
dt 



A 



iV(x-£,y)(l-^%^)-7V(x,y) 



dx 



E 



e{E) - + D AK 



where O is the Heavyside step function and 
E = \Vh\ - tan{Or). 



(7) 



(8) 



The new model is called the Saltation-Creep- Avalanche 
(SCA) model. It will be investigated in the present paper 
and numerical results will be compared with respect to 
NO landscapes. One should note that we have numeri- 
cally solved Eq.(7) using the same flux law (5) as in the 
original NO model. 



One observes that : (i) the growth of the maximum 
height hmax is linear and (ii) there is no saturation. This 
means that ripples become inflnitely high! This unreal- 
istic behavior comes from the right hand side of Eq. (4) . 
Indeed, when the flrst term becomes larger, the sand de- 
posit becomes more important than the relaxation. As 
time goes on, a net growth of the proflle is thus recorded. 
This situation will remain unchanged since the ratio be- 
tween deposit and relaxation is fixed once by the values 
of 6, e and D. Consequently, hmax grows ad infinitum. 
One should note that this pathology implies the possi- 
bility of making cubes of sand since nothing prohibits 
vertical walls. In order to get a more realistic situation, 
an additional term should be included in the NO model. 
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FIG. 2. Avalanches lead to a profile stabilisation: (left) 
the local slope tan(S) is larger than the critical value tan(0j.), 
the surface is unstable, (right) The avalanche brings the 
local slope to a new stable value tan{9r). Note that 
E = IVftj - tan(6l^). 
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III. RIPPLE FORMATION 

In this section, we investigate the influence of the four 
parameters S, e, D and 0,. on the formation of the rip- 
ples in the SCA model. All the simulations presented 
herebelow are made on 150 x 150 square lattices. The 
evolution of the landscape is typically stopped after 10^ 
iterations. The default parameter values are S = 12.5, 
£ = 0.3, D = 0.4 and Or = 37°. 

A. Flux parameter S 

The saltation process obeys Eq.(5) in which the pa- 
rameter 6 controls the flux Qsait- A small arbitrary 
value of 5 implies that small amounts of sand are dis- 
placed along small saltation lengths, while larger values 
of S lead to a more important transport of sand over 
long distances. Experimentally [1], one observes that the 
ripple wavelength (the distance between two consecutive 
crests) grows with the saltation length £. 

We have performed simulations varying S from 1 to 
100. For very small values of the parameter {6 < 7), no 
ripple appears and the surface remains irregular. This 
threshold is not encountered in the NO model and is 
certainly due to the angle of repose. For small values of 
the flux parameter, not enough sand is displaced and the 
surface is continuously smoothed by avalanches. 

Above this threshold S ~ 7, a set of several ripples 
appears. The density of ripples depends on 6. The larger 
5 is, the smaller is the density of ripples. For example, 
after 5 10® iterations, the density of ripples is 0.04 for S = 
10, while there is only a single ripple (density of 0.006) 
if 5 = 85. The SCA model is thus in good agreement 
with observations: the wavelength of the ripple pattern 
grows with the saltation length i. In the case of the NO 
model, the density of ripples grows non-linearly with S in 
opposition with experimental observations. 
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FIG. 3. The normalized saltation flux Q salt/ 5 of Eq.(5) as 
a function of the spatial derivative Vh of the height of sand 
h. Different values of the parameter e are illustrated. 



Practically, we have considered values of s ranging 
from to 2. Our observations are: (i) no ripple appears 
if £ = 0. (ii) for small values of e (typically < £ < 0.2), 
ondulations appear in the landscape, (in) As e becomes 
larger, the number of ripples present in the system grows 
(from 5 at £ = 0.5 to 11 for e = 1.8). Moreover, rip- 
ples are more asymmetric (windward slopes smaller than 
downstream slopes) for large values of e. We will study 
this asymmetry in section 4. 

In the case of the NO model we observe that: (i) if 
e = very small ripples appear, (ii) when the parameter 
reaches a threshold e ~ 0.2 larger ripples appear. Note 
that the wavelength of the ripples does not depend on s 
above that threshold. The fact that ripples appear for 
£ = in the NO model and do not appear within the 
SCA model is due to the existence of an angle of repose. 
Indeed, avalanches tend to smooth the landscape. 

In summary, the flux Qsalt should be asymmetric in 
order to create ripples. For positive slopes, a minimum 
value of this flux is required. 



B. Asymmetry parameter e 

In the SCA model, we have assumed a non-zero value 
for the flux of sand Qsalt displaced on the stoss slopes 
(which are upstream) of the ripples (see Eq.(5)). This 
hypothesis is supported by numerical results [3]. If the 
parameter e is non zero, the saltation flux becomes asym- 
metric (see Figure 3). Sand is always taken from the 
stoss slope, while Qsalt tends to zero on the lee slope. 
This comes from the fact that sand on the lee slope is 
confined in wind vortices. 



C. Reptation coefficient D 

The parameter D plays the role of a diffusion coefficient 
for smoothing the irregularities along the surface. We 
have considered values of D ranging from to 1. Once 
again, no ripple appears for very small values of that 
parameter. Ripple appearance is allowed in the interval 
0.1 < D < 0.5. If £> < 0.1 the surface remains irregular. 
There is nearly no interaction between the cells because 
of the small diffusion coeflicient. Note that the same 
behavior is observed in the NO model. As D becomes 
larger, the ripple wavelength increases. Finally, for D > 
0.5 the diffusion coeflacent becomes too large and the 
surface remains smooth and flat. For the NO model, 
no ripple appears when D > 1.8, a larger value than in 
the case of the SCA model. 
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D. Angle of repose Or 

Parameters D and Or are intrinsic parameters of the 
granular material, while 5 and e shoiild be related to the 
wind. At each iteration, the reptation acts and moves 
grains. On the other hand, an avalanche occurs if and 
only if the local slope of the surface exceeds the angle of 
repose. The role played by both D and Or parameters is 
thus different: D is a diffusion coefficient, which controls 
the flow of sand between neighboring cells, while the an- 
gle of repose Or governs the local slope of the granular 
surface. 

We considered typical values of Or between 20° and 
40°. For ripples of the same wavelength, we have noted 
that small angles of repose imply small ripple heights. 
The angle of repose is thus responsible of the value of the 
ratio hmax/\ where A is the ripple wavelength. Indeed, 
a slope larger than tan(6'r) is reduced to a value tan(6'r) 
by avalanches. We have also noted the appearance of 
ripples for all tested values of Or- 



IV. RIPPLE DYNAMICS 

In this section, we investigate the dynamics of the SCA 
model. Three different aspects will be considered: the 
ripple height, the ripple shape and the formation of kinks. 
Again, the default parameters values are 5 = 12.5, e = 
0.3, D = 0.4 and 61^ = 37°. 



A. Ripple amplitude hmax 

Figure 4 presents the temporal evolution of the max- 
imum ripple height h„iax in the case of the SCA model. 
One should note the saturation of hmax for long times. 
This figure should be compared with Figure 1. 
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FIG. 4. Evolution of the maximum height hmax of the land- 
scape as a function of time t in the case of the SCA model. 
The data have boon fitted using Eq.(9). The fit parameters 
are a = 230.8 ±0.6, b = 377.3 ±3 and c = 1.33 lO'^iO.Ol lO'^. 
The simulation parameters are the sames as those of Figure 
1. The angle of repose is Or = 37°. When all ripples have 
merged, the system contains only one ripple for which the 
height is limited by the angle of repose. 



In Figure 4, many gaps can be observed in the early 
stages of evolution. Those gaps come from the absorp- 
tion of small ripples by the largest ones. Since a small 
ripple travels faster than a large one [1] , the large ripple 
tends to absorb the small. Figure 5 presents the evo- 
lution of a part of the lattice for different stages of the 
simulation. For each picture, a transverse view of the 
landscape is shown. One should note that the merge 
occurs as follows: first, the small ripple climbs on the 
larger one. The arrival of the small ripple on the crest 
causes an avalanche. A part of the small ripple falls in 
the avalanche, while the remainder feeds the larger one. 
It results in a fast growth of the largest ripple even it 
moves slowly. 
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FIG. 5. The merge of two ripples. The wind blows from 
left to right. The profile of tho landscape is shown below 
each picture, (top row) the smallest ripple travels faster than 
the largest one. (second row) The small ripple climbs on the 
large one. (third row) An avalanche is initiated by the arrival 
of the small ripple on the crest, (bottom row) The small ripple 
'feeds' the large one. Note that the profiles have been rescaled 
in order to be compared. 
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Asymptotically, only one ripple occupies the whole lat- 
tice, and the time needed for a global merge becomes 
larger as S decreases. If S is large, the initial density of 
ripples is small. In this case, the global merge is obtained 
after a small number of collapses. In opposition, a large 
number of merges is necessary for small values of S. If one 
assumes that ripple collapses is a size independent pro- 
cess, one can understand that the global merge is faster 
when (5^1. 

The main characteristics of hmax{t} '^^^ (i) S^ps in the 
early stages of ripple formation, and (ii) & saturation for 
long times. Looking for details in Figure 4, one could 
see that gaps are allways followed by the same kind of 
growth. Indeed, the maximum ripple height evolves ac- 
cording to an exponential growth law 



b exp 



(9) 



where a, b and c are fitting parameters. This law is shown 
in Figure 4. One should see that the fit has been per- 
formed after the primary gaps, e.g. for t > 10^ iterations. 

In summary, we have seen that, contrary to the NO 
model, the SCA model naturally leads to a non-linear 
evolution of hmax- Note that this kind of behavior has 
been experimentally observed in [10,11]. Moreover, the 
SCA model predicts a saturation of hmax after a finite 
time, as expected and observed. The saturation value 
depends essentially on 0r. The larger is 0r, the higher 
are the ripples. 



B. Ripple aspect ratio a 

Under a non-oscillating wind, a ripple has generally an 
asymmetric shape [1,8]. The stoss slope is indeed smaller 
than the lee slope. In Figure 6, the profile of the ripples 
is shown for both NO and SCA models. Ripples are more 
symmetric in the NO model than in the SCA model. 



In order to measure the aspect ratio cr of a ripple, we 
have measured the ratio of the horizontal projections of 
both sides of the ripples. Let us call Xg the projection 
of the stoss slope on the horizontal x axis, and xi the 
projection of the lee slope. The apect ratio of a ripple is 
defined as 



Xg 



(10) 



The measure of a is averaged over all the ripples of the 
landscape. If cr = 1 the ripples are symmetric, while 
cr ^ 1 implies an asymmetry. 

In Figure 7, the temporal evolution of a for both NO 
and SCA models are reported. Typical values of cr range 
between 1 and 5 in the SCA model. In the NO model, a 
remains close to 1. 

In the case of the SCA model, one should note that 
cr grows linearly with time, and then remains close to a 
saturation value. The appearence of a saturation seems 
coherent since only one ripple occupies the whole lattice 
after a finite time. The angle of repose limits the height 
of that ripple and its morphology. 
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FIG. 7. Temporal evolution of the aspect ratio a. Both 
NO and SCA models are illustrated. The parameter values 
are 5 = 12, e = 0.3, D = 0.4 and Or = 37°. 
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FIG. 6. Height h of sand (in arbitrary units) as a function 
of the spatial coordinate x. Both NO and SCA models are 
shown. The curves have been rescaled in order to be com- 
pared. 



In order to compare the four parameter dependencies 
on (7, we have considered a reference situation which cor- 
responds to the default parameters. Varying the value of 
a parameter while keeping constant the others allows us 
to define the infiuence of that parameter on the temporal 
evolution of the aspect ratio. Herebelow, we describe the 
results obtained for each of the four parameters S, e, D 
and 6r- One should note that in all cases, the temporal 
evolution of a begins with a quasi-linear law, followed by 
a saturation. 

The aspect ratio increases with the parameter S. Since 
the ripple wavelength directly depends on that result 
means that small ripples are less asymmetric than large 
ones. This result is in good agreement with observations 
[1,8]. 

Figure 8 presents the temporal evolution of the aspect 
ratio C7 for different values of s. One can see that small 
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ripples (for small e values) are less asymmetric. Indeed, 
(7 is larger for e = 1 than for e = 0.3. We have seen in 
Section III.B. that the ripple v^favelength becomes larger 
as e increases. One should also note that a saturates 
more rapidly when e is small. This is consistent with the 
fact that the time needed for a global merge is shorter in 
the case of large ripples. 

As observed for 5 and e, the aspect ratio follows any 
variation of D. This results is in agreement with those 
observed previously in this section. Indeed, we have seen 
in Section III.D. that the ripple wavelength increases as 
D becomes larger. 

We have taken realistic values of Or between 30° and 
40°. The aspect ratio a is found to be independant of 
the value of Or- This result is outside the scope of this 
paper and will be studied in the future. 
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FIG. 8. Temporal evolution of the aspect ratio a in the 
case of the SCA model. Different values of e are illustrated. 



C. Kink dynamics 

Although the wind blows in a single direction, the re- 
sulting ripples are not strictly perpendicular to the wind 
direction. There are defects. Two different classes of de- 
fects can be found: kinks which are fusion points of two 
crests, and antikinks which are endings of crests. In Ref. 
[7] , wc have studied in the NO model how defects modify 
a landscape morphology when a brutal change occurs in 
the wind direction. We have also shown in Ref. [7] that 
under a constant wind direction, the density of defects 
decreases as time goes on. In the present section, we will 
focus on the processes leading to this decrease in the SCA 
model. 

After extensive simulations, we conclude that some 
kinks can propagate during very long times, while others 
disappear quite rapidly. 

The disappearance of a defect mainly occurs when a 
ripple is splitted into two inequitably parts (see an ex- 
ample in Figure 9). The small branch moves faster than 
the larger one and climbs on it following the merging dy- 
namics described in section IV. A. Reaching the crest, the 



small branch creates an avalanche. The sand rolling on 
the lee slope may be captured. If the amount of sand 
escaped from this avalanche is small, the defect disap- 
pears and the small branch vanishes. After this process, 
the number of defects is reduced. If all defects show this 
kind of behavior the defect density will rapidly tend to 
zero. However, this is not the case [7]. 

The propagation of a kink takes place when a sufS- 
cient amount of sand escapes from the avalanche. This 
is the case when the branches are quite long. In the em- 
phasized part of Figure 10, one can see the ending of 
a small ripple. The extremity of the ripple moves faster 
than the main part. The merge (second image) is accom- 
panied by an avalanche. As seen in the previous para- 
graph, sand is expelled out of the avalanche. But here, 
the amount of sand is sufficient to create a new ripple. 
The merging/propagation process is thus repeated on the 
next ripple. 

Let us consider a branch of length ib. The merge of 
this branch with a larger ripple causes the decrease of the 
branch length since a part of the branch has been owned 
by the ripple. When a second meeting occurs, ib is once 
again reduced. As a consequence, the initial length of 
the branch is a relevant parameter in order to predict 
the life-time of a defect. For example, a kink of initial 
length £b = 65 can typically propagate during 9.5 x 10^ 
iterations, while a kink of initial length £b = 23 has a life 
expectancy of 2.5 x 10^ iterations. One should note that 
those life-times mainly depend on the diffusion coefficient 
D and the avalanche process. 

In Ref. [7] , we have found a non zero asymptotic value 
for the defect density. The existence of propagating kinks 
gives a picture of this behavior. The relevant ingredient 
for kink disappearance being the length of the branch. 




FIG. 9. The disappearance of a kink, (top row) The small- 
est branch of the ripple climbs on the main part, (bottom 
row) the avalanche feeds the ripple, and the branch vanishes. 
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FIG. 10. The propagation of a kink. The defect is expelled 
downstream. 



V. SUMMARY 

We have introduced new ingredients in the Nishimori- 
Ouchi (NO) model: the existence of an angle of repose 
Or and subsequent avalanches when the local slope 9 is 
larger than Or- The new model reproduces realistic fea- 
tures of aeolian ripples such as a non-linear evolution of 
the ripple amplitude hmax- The origin of this behavior 
has been explained by the merge of ripples traveling at 
different speeds. Increasing the saltation length, we have 
observed a grows of the ripple wavelength, in agreement 
with observations. In the case of a constant wind orien- 
tation, natural ripples are asymmetric. This feature has 
also been reproduced by the new model. 

Studying the kink dynamics, we have shown the coex- 
istence of two kinds of kink evolution: propagation and 
disappearance. The defect initial length has been pro- 
posed as a relevant parameter in order to predict the 
life-time of that defect. From this result, the existence of 
a non-zero asymptotic value of the kink density has been 
explained. 
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